clear, clc 
global w v 
% 设置 水分生态位和温度
w = zeros(3, 1); w(1) = 1.488; w(2) = 1.75; w(3) = 3.19; 
v = zeros(3, 1); v(1) = 17.82; v(2) = 15.32; v(3) = 19.9; 

% 设置 相对水分生态位和温度
w = w./sum(w); v = v./sum(v); w = w'; v = v'; 
wa = -0.2; T = 26; 
xmax1=2000;
xmax2=2000;
xmax3=2000;
x10=1;
x20=1;
x30=1;
r1=0.8;

%[t, x] = ode45(@(t,x)func(t,x,wa,T), [0 150], [100 200 200]);
% 单物种无环境变化图no environmental change for single species
t=60
[t, x] = ode45(@(T,x)func1(t,x,wa,T), [0 150], [100 200 200]);
noE4Single = figure('Name','no environmental change for single species','NumberTitle','off');
plot(t,x)
title('no environmental change for single species')
xlabel('time(day)')
ylabel('n')
saveas(noE4Single,"../../figures/noEChanged4SingleSpecies.png")

% 单物种有环境湿度变化图moisture change for single species
t=60
[t, x] = ode45(@(T,x)func22(t,x,wa,T), [0 150], [100 200 200]);
mositureChanged4SingleSpecies = figure('Name','moisture change for single species','NumberTitle','off');
plot(t,x)
title('moisture change for single species')
xlabel('time(day)')
ylabel('n')
saveas(noE4Single,"../../figures/mositureChanged4SingleSpecies.png")

% 单物种有环境温度度变化图temp change for single species
t=60
[t, x] = ode45(@(T,x)func23(t,x,wa,T), [0 150], [100 200 200]);
tempChanged4SingleSpecies = figure('Name','temp change for single species','NumberTitle','off');
plot(t,x)
%标题
title('temp change for single species')
% 标签
xlabel('time(day)')
ylabel('n')
saveas(noE4Single,"../../figures/tempChanged4SingleSpecies.png")

% 多物种有环境变化图environmental change for single species
t=60
[t, x] = ode45(@(T,x)func32(t,x,wa,T), [0 150], [100 200 200]);
tempChanged4SingleSpecies = figure('Name','E change for multi species','NumberTitle','off');
plot(t,x)
%标题
title('E change for multi species')
% 标签
xlabel('time(day)')
ylabel('n')
saveas(noE4Single,"../../figures/Changed4multiSpecies.png")



















[t, x] = ode45(@(t,x)func23(t,x,wa,T), [0 150], [100 200 200]);

[t, x] = ode45(@(t,x)func32(t,x,wa,T), [0 150], [100 200 200]);

x1=xmax1./(1+((xmax1-x10)/x10)*exp(-r1*t));
plot(t,x1,'r-')

hold on

r2=0.4;
x2=xmax2./(1+((xmax2-x20)/x20)*exp(-r2*t));
plot(t,x2,'g-')

hold on

r3=0.2;
x3=xmax2./(1+((xmax2-x30)/x30)*exp(-r3*t));
plot(t,x3,'b-')

%x=diff(x1);
%td=diff(t);
d1=2.518892+0.2980319.*(r1-x1./xmax1).*(r1-x1./xmax1)+3.059665.*exp(0.877)
plot(d1)

hold on
%x=diff(x2);
%td=diff(t);
d2=2.518892+0.2980319.*(r2-x2./xmax2).*(r2-x2./xmax2)+3.059665.*exp(0.696)
plot(d2)
hold on

%x=diff(x3);
%td=diff(t);
d3=2.518892+0.2980319.*(r3-x3./xmax3).*(r3-x3./xmax3)+3.059665.*exp(0.315)
plot(d3)

ht=r1-r1.*x1/xmax1;
d1=2.518892+0.2980319*ht*ht%+3.059665*exp()